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ABSTRACT 


An  analog  system  subject  to  the  Poisson  Shock  is  modeled 
using  past  performance  data.  Failure  Dynamics  of  the  system 
is  estimated  by  curve  fitting  techniques.  Algorithms  for 
fault  prediction  in  an  on-line  maintenance  process  are  describ 
ed.  Several  sequential  refinement  schemes  are  introduced  to 
improve  fault  prediction.  Some  formulas  and  properties  of 
system's  statistics  have  been  developed.  A decision  rule  is 
introduced  which  is  based  on  the  criteria  of  simultaneously 
maximizing  lifetime  and  minimizing  the  cost  of  on-line  fail- 
ures. Poisson  Shock  generator  is  implemented  by  computer  for 
simulation  of  the  on-line  maintenance  process.  The  computer 
simulations  of  a perfect,  no  measurement  errors  and  identical 
drifting  parameters,  system  are  presented.  The  simulations 
of  an  imperfect  system  are  studied  by  adding  a noise  to  the 
system  performance  data. 
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INTRODUCTION 

Fault  analysis  processes  have  been  and  will  continue 
to  be  very  significant  factors  in  the  safety  and  reliability 
of  electrical  systems.  This  is  especially  true  due  to  the 
following  facts:  a rapid  advancement  in  the  complexity  and 
size  of  modern  systems;  increased  availability  and  capabili- 
ties of  computers;  and  rapidly  changing  technologies  in  in- 
tegrated circuit  fabrication.  Due  to  this,  fault  analysis 
has  become  much  more  than  an  academic  research  topic.  Fault 
analysis  is  applicable  in  an  industrial  environment  to  mini- 
mize cost,  extend  the  lifetime  of  the  overall  system,  con- 
trol maintenance  schedules,  and  effectively  plan  manpower 
needs . 

The  techniques  of  fault  analysis  (detection,  location 
and  prediction)  have  been  developed  independently  in  two 
different  areas,  digital  and  analog  systems.  In  the  last 
decade,  considerable  fault  location  research  has  been  per- 
formed in  the  digital  system  area.  Preparata,  Metze  and 
Chien  (1)  introduced  an  elegant  graphical  model  and  defined 
the  properties  of  t-diagnosable  systems.  This  model  proved 
to  be  very  powerful  with  the  development  of  characterization 
theorems  for  the  system  by  Hakimi  and  Amin  (2,  3).  In  1973, 
Russell  and  Kime  (4)  extended  the  Preparata -Met ze-Chien 
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model  into  a more  general  case.  However,  such  a generali- 
zation increases  the  complexity  of  the  system. 

Fault  detection  research  for  digital  systems  was 
introducted  by  Patterson  and  Metze  (6)  in  1972  and  by 
Pradhan  and  Reddy  (7)  in  1973.  Several  different  methods 
were  developed  for  solving  the  problem  of  real-time  fault 
detection  in  an  asynchronous  sequential  circuit.  After 
1974,  Maki  and  Sawin  III  (8,  9,  10)  published  a series  of 
papers  discussing  the  fault  detection  and  fault-tolerant 
design  for  both  synchronous  and  asynchronous  sequential 
circuits. 

The  use  of  graphical  models  for  dealing  with  analog 
circuits  was  first  suggested  by  Kirchhoff  in  1847  (11) . 
Using  Kirchhoff' s model,  Berkowitz  (12)  developed  his  con- 
cept of  solvability  for  fault  diagnosis  in  1960.  From  the 
necessary  conditions  for  solvability,  a new  concept.  Key 
Subgraph  (13),  was  introduced.  Gayer  (14)  applied  the  Key 
Subgraph  concept  to  fault  isolation  in  simple  linear  solid 
state  amplifiers. 

Another  aspect  in  the  fault  analysis  process  for  an 
analog  device  is  fault  prediction.  To  accurately  predict 
a fault,  a device  must  be  tested  at  periodic  maintenance 
intervals.  If  the  device  fails  or  does  not  operate  cor- 
rectly, it  is  replaced  immediately.  The  device  may  be 
assumed  good  if  its  characteristics  are  normal.  However, 
if  the  characteristics  are  slightly  out  of  tolerance,  but 
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the  device  still  operates  correctly,  one  can  attempt  to 
predict  if  the  device  will  fail  before  the  next  scheduled 
maintenance  interval.  If  device  failure  is  predicted,  it 
can  be  replaced  before  failure  occurs  as  part  of  planned 
preventive  maintenance. 

With  the  advent  of  the  low-cost  microprocessor,  on- 
line fault  prediction  is  possible  and  practical.  A curve 
fitting  algorithm  for  on-line  fault  prediction  was  first 
introduced  by  Saeks,  Liberty  and  Tung  (15,  16,  17)  in  1975. 

It  was  assumed  that  prior  life-time  statistics  for  the 
system  under  test  were  known.  Also,  performance  data  of 
the  system  at  each  maintenance  interval  were  collected. 

The  application  of  these  data  to  a second  order  polynomial 
equation  resulted  in  an  estimation  of  the  time  at  which 
the  component  under  test  would  exceed  tolerance  limits. 

Based  on  a criterion  of  simultaneously  minimizing  on-line 
failures  and  maximizing  component  lifetime,  a decision  as 
to  whether  or  not  the  component  should  be  replaced  is  made 
at  each  maintenance  interval. 

The  disadvantages  of  this  curve  fitting  algorithm  are: 
the  application  is  limited  to  failures  due  to  permanent 
over stress,  the  second  order  polynomial  is  too  simple  to 
describe  the  performance  of  the  component,  and  the  prior 
lifetime  statistics  for  the  component  are  often  not  available. 

Another  area  where  an  extensive  research  effort  is 
being  made  is  shock  models  and  wear  processes.  Esary,  Mar- 
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shall  and  Proschan  (21,  22,  23)  introduced  a shock  model 
for  the  life  distribution  of  a component  subjected  to  a 
sequence  of  shocks  randomly  occurring  in  the  time  according 
to  a homogeneous  Poisson  process.  They  also  considered 
the  related  shock  models  in  which  each  shock  caused  a 
random  amount  of  damage  and  failure  occurred  when  the 
accumulated  damage  exceeded  a specified  threshold.  This 
failure  model  is  well  known  in  modern  reliability  theory. 

Employing  the  Poisson-Shock  model,  another  curve 
fitting  fault  prediction  algorithm  which  will  overcome  the 
disadvantages  of  the  Saeks-Liberty-Tung  algorithm  will  be 
discussed  in  the  succeding  chapters.  In  Chapter  2,  an  ana- 
log system  subjected  to  the  Poisson  Shock  is  modeled  using 
past  performance  data  of  the  analog  system.  Algorithms  for 
fault  prediction  in  an  on-line  maintenance  process  are  de- 
scribed. Several  sequential  refinement  schemes  have  been 
developed  to  improve  fault  prediction.  In  Chapter  3,  a 
decision  rule  is  developed  which  is  based  on  the  criterion 
of  simultaneously  maximizing  lifetime  and  minimizing  the 
cost  of  on-line  failures.  Simulations  of  the  proposed  on- 
line maintenance  process  are  presented  in  Chapter  4. 
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CHAPTER  2 
SYSTEM  MODELING 


2.1  Failure  Dynamics 

The  performance  of  an  analog  device  subject  to  a 
series  of  discrete  shocks  (switching  process,  improper 
operation,  etc...)  may  drift  due  to  the  shock  damage.  Let 
C (N)  represent  values  of  a particular  component  parameter, 
where  the  component  time,  N,  denotes  the  number  of  shocks 
the  component  has  received.  In  order  to  normalize  the 
value  of  parameter  C,  a 1 or  0 is  assigned  to  a "brand  new" 
or  "fail"  stage  respectively. 

The  normalized  value  of  the  parameter  should  have 
the  following  properties: 

(I)  C (0)  = 1 

(II)  C (°°)  = 0 (2.1-1) 

(III)  If  L is  the  smallest 
integer  that  C(L)  =*  0,  then 
C(N)  = 0 V N - L . 

It  is  assumed  that  drifting  parameters  can  be  well  described 
by  a difference  equation  of  the  following  form: 

C(N  + 1)  = C (N)  - f (N)  N < L (2.1-2) 

where  f (N)  is  the  particular  failure  dynamics. 
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Due  to  the  simplicity  of  the  calculation,  it  is  practi- 
cal to  assume  that  f(N)  has  a polynomial  form  with  order  h. 
That  is,  eq.  (2.1-2)  can  be  rewritten  as 


C(N)  - C (N+l)  * aQ  + axN  +....+  ahN 


h i 

,z  aiN 

1*0 


(2.1-3) 


where  N < L and  the  {a^}  are  constants  and  L is  the  smallest 
integer  satisfying  the  condition 


L-l  h 

Z Z a . j1  > 1 

j*0  i-0  1 ” 


(2.1-4) 


Solving  the  difference  eq.  (2.1-3)  with  the  boundary  condi- 
tions (2.1-1),  one  obtains 


N-l  h 


C (N)  * 1 - Z 


Z a,  j‘ 
j*0  i-0  1 


C(N) 


if  N < L 


if  N > L 


(2.1-5) 


L is  termed  the  "life-time"  of  the  component,  that  is,  the 
number  of  shocks  necessary  to  cause  the  component  to  fail. 

Consider  a simple  example  where  f(N)  is  taken  to  be 
a first  order  polynomial;  that  is, 


f (N)  - C(N)  - C (N+l)  * aQ  + ajN 


(2.1-6) 

Prom  equation  (2.1-6)  and  boundary  conditions  (2.1-1), 
C(N)  can  be  expressed  as 
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C(N)  - 1 - aQN  - N-^—  ax  if  N < L 
C(N)  - 0 if  N > L 


(2.1-7) 


Then  the  life-time  of  this  component  is  the  smallest  integer 
satisfying  the  equation 


1 - aQL  - aL  > 0 (2.1-8) 


That  is. 


L - 


> - a^2  + 8a^  - (2aQ  - a^ 


(2.1-9) 


2a, 


Table  (2.1-1)  gives  the  life-time  of  the  component, 
satisfying  eq.  (2.1-8),  with  different  values  of  aQ  and  a^. 
Fig.  (2-1-1)  shows  the  life-time  curve  as  a function  of  aQ 
with  different  a^'s. 


Component's  Lifetime  with  different  values  of  art  and 
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2.2  Poisson  Shock 

Shocks,  received  by  the  component,  may  occur  randomly 
in  time  according  to  some  statistical  process.  The  Poisson 
process  is  a well-founded  shock  model  in  modern  reliability 
theory.  That  is,  the  probability  of  N shocks  occuring  in 
the  time  interval  t is: 


Pt(N) 


e-kt  (kt)N 
N! 


N = 0,  1,  2 


(2.2-1) 


Where  k is  a given  constant  representing  the  average  number 
of  shocks  per  unit  time.  Therefore,  (kt)  is  the  average 
number  of  shocks  in  the  time  interval  t. 

Let  F ( t)  be  denoted  as  the  probability  a component 
subjected  to  a Poisson  shock  with  constant  k,  will  not  fail 
by  time  t.  Then, 

L-l 

F ( t)  = E p.  (n) 
n=0 


K1  _-kt  (kt) n 
■ l e 

n=0  n 1 

Table  (2.2-1)  gives  the  solution  of  F(t) 


(2.2-2) 

for  different 


values  of  L and  (kt) . Fig.  (2-2-1)  shows  F(t)  as  a function 
of  (kt)  with  different  L's. 


Solutions  of  F(t) for  different  values  of  L and  (kt) 
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2.3  Estimation  of  the  Failure  Dynamics 

In  a periodic  maintenance  system,  the  performance 
data  of  a component  is  taken  at  each  maintenance  interval 
nT.  Thia  is,  (C^,  ....,  C^)  is  the  performance  data 

taken  at  maintenance  times  (T,  2T  ....,  gT)  . The  problem 
can  be  states  as: 

"Given  performance  data  (C^  C2  ....,  Cg)  , solve 
for  the  constants  (aQ,  a^  ...f  a^)  of  the  failure  dynamics 

in  eq.  (2.1-3) . " 

Since  it  is  assumed  that  the  system  is  subjected  to 
Poisson  Shock  with  constant  k,  that  is  the  average  number 
of  shocks  in  each  time  interval  T is  kT.  Therefore,  eq. 
(2.1-5)  can  be  rewritten  as: 


mkT-1  0 mkT-1  ^ mkT-1  h 


E a0^  + 


E “1 


a.j  + + 


Z h 


■ 1 - C 


where  m=  1,  2,  3,  ....,  g 
or  in  matrix  form: 


(2.3-1) 


kT-1 

kT-1 

kT-1 

Z 

j° 

E 

j1 

....  E 

i_i. 

11 

0 

j-0 

j-0 

2 kT-1 

A 

2 kT-1 

1 

2kT-l 

Z 

j° 

E 

j1 

....  E 

j-o 

j-0 

j-0 

gkT-1 

0 

gkT-1 

1 

gkT-1 

E 

E 

j1 

....  E 

j-0 

j-0 

j-0 

l '] 

(2.3-2) 


J * A = Z 


(2.3-3) 


where 


kT-1 

E 

3° 

kT-1 

E 

31  • • • 

kT-1 
. E 

j-o 

j-o 

j=o 

2 kT-1 

2kT-l 

1 

2kT-l 

E 

j° 

E 

j • • • 

E 

j-o 

j-o 

j-o 

• • • 

gkT-1 

• • 

o 

• • • 

gkT-1 

• • • • • 

1 

• • • • 

gkT-1 

E 

3 

E 

j • • • 

E 

j-o 

o 

H 

•n 

j-o 

Z » l-c 


A = a. 


The  solution  A cam  be  discussed  in  two  different  cases: 
(i)  If  J is  nonsingular,  then 

A - J-1  * Z (2.3-4) 

(ii)  If  J is  singular  and  full  column  rank,  then 


A - (JTJ)_1  JTZ 


(2.3-5) 
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2.4  Sequential  Refinement  Schemes 

To  approximate  the  failure  dynamics  by  a polynomial, 
the  proper  choice  of  the  order  h is,  in  general,  quite 
difficult  and  depends  upon  physical  considerations  and  the 
engineering  experience.  Once  h is  preselected,  employing 
the  technique  developed  in  the  previous  section,  the  poly- 
nomial coefficients  to  best  approximate  the  failure  dynamics 
cam  be  calculated.  The  accuracy  of  the  failure  dynaunics 
depends  greatly  on  the  choice  of  the  order  and  can  be  im- 
proved by  considering  more  measurements.  To  find  a new  set 
of  coefficients  for  a different  combination  of  h and  g, 
the  entire  calculation  procedure  must  be  repeated  from  the 
very  beginning.  This  repitition  is  impractical  in  the  on- 
line maintenamce  system.  Therefore,  sequential  refinement 
schemes  for  obtaining  new  sets  of  coefficients  without 
having  to  repeat  the  entire  calculation  will  be  developed 
for  the  case  in  which  g and  h are  allowed  to  vary.  It 
is  assumed  that  the  coefficients  matrix  A,  found  from  the 
last  section  for  an  h order  polynomial  with  g measure- 
ments, are  solutions  for  eq.  (2.3-2).  Three  typical  cases 
will  be  considered  in  the  succeeding  paragraphs. 
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Case  1: 


Suppose  that  an  additional  measurement  is  taken 
at  the  next  scheduled  maintenance  interval  (g+l)T,  then  eq. 
(2.3-2)  becomes: 


kT-1 

kT-1 

kT-1 

E 

3° 

E 

j1 

. . . E 

jh 

an 

i-C 

j-o 

j-o 

j-o 

u 

1 

2kT-l 

E 

j° 

2kT-l 

E 

.1 

3 

2 kT-1 
. . . E 

ih 

al 

1-C2 

j-o 

j-o 

j“0 

JL 

• 

• 

gkT-1 

o 

gkT-1 

1 

gkT-1 

h 

• 

• 

E 

j 

E 

j 

. . . E 

3 

1“C_ 

j*0 

j-o 

j— o 

• 

9 

ah 

(g+1) kT- 

i0 

(g+1) kT- 

1 1 

(g+1) kT- 

1 h 

n 

E 

j 

E 

j • 

. . E 

3 

1-cg+i 

j-o 

j-o 

j-o 

(2.4-1) 


or 


(2.4-2) 


where  J and  Z were  defined  in  eq.  (2.3-3) 


(g+1) kT-1  n 

(g+1)  kT- 

■1  . (g+1) kT- 

E j° 

E 

j\  . . E 

j“0 

j-o 

j-o 

A 

A * 

[a0  al  * * ‘ ah^ 

A 

Z 

1-C  .. 

g+1 

mm 
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The  matrix  A cam  be  solved  as : 


rp  .1  m 

A = (J  J + J J)  1 ( j Z + J Z)  (2.4-3) 
To  avoid  the  calculation  of  a new  matrix  inverse  as  indica- 


ted in  eq.  (2.4-3)  a sequential  algorithm  (24,  25)  may  be 
used  to  update  A.  That  is,  assume 


T -1 
P - (J  J)  1 


is  given,  then  from  eq.  (2.3-5) 


A » PJ  Z 


and , 


m A«pA  .1  Am  ^ Afp  .1  A 

(J  J + J J)  * P - PJ  (JPJ  + I)  ■LJP  (2.4-4) 


Therefore, 


A - P - pj 


a *¥•  A A a 1 A IP  A IP  A 

^(jpj  ♦ :)  AJP(J  z + j z) 


(2.4-5) 


A»p  A aIP  • 1 A A rn  A AlP  A a fP  M 1 A 

a + pj1  { r-(jpjA+i)  1jpj1}  z - pj  (jpj  +i)  xja 


using  the  identity 


~ - I O’  • V * •'•m  _ 1 

1 - (JPJ  ♦ I)  ^JPJ1  * (JPJ1  + I)  1 


(2.4-6) 


Eq.  (2.4-5)  becomes 

A A IP  * A IP  a,  1 A A 

A - A + PJ*(JPJ  ♦ I)  (S-JA)  . 


(2.4-7) 


Notice  that  the  manipulation  of  the  matrix  (JPJT  + I)~^ 


is  a simple  scalar  inversion. 


t 


1 
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Case  2: 

Suppose  an  h order  polynomial  has  been  used  and  it 
is  then  determined  that  the  error  is  too  large.  It  would  be 
more  desireable  to  use  a higher  order,  H,  polynomial  without 
having  to  repeat  the  entire  calculation.  Then,  eq.  (2.3-2) 
becomes : 


' kT-1 

E 

j° 

kT-1 

E 

.1 

3 • 

kT-1 
. E 

.h 

3 

kT-1 
. . E 

• H 

3 

j-o 

j = 0 

j=0 

LJ. 

II 

o 

2 kT-1 

E 

.o 

3 

2kT-l 

E 

j1  . 

2kT-l 
. 1 

.h 

3 

2kT-l 
. . E 

.H 

3 

j-0 

j-o 

j-0 

j-0 

gkT-l 

E 

3° 

gkT-l 

E 

j1 

gkT-l 

E 

.h 
3 1 

gkT-l 

E 

.H 

3 

j-o 

j-o 

j-0 

j-0 

- 

■ ■ 

■ 

a 

o 

L-Cl 

al 

• 

• 

s 

1“C2 

• 

ah 

• 

• 

> 

1-C 

g 

or 


A 

J 

J] 

/N 

B 

(2.4-8) 


(2.4-9) 


where  J and  Z were  defined  in  eq.  (2.3-3), 


kT-1 

E j 
j-0 


h+1 


kT-1 
E 3 

j-o 


H 


j=o 


gkT-l 

£ jH 

• J-0 
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In  order  to  improve  the  computational  efficiency,  an 
algorithm  {24,  25)  may  be  used  to  compute  the  matrix  inverse 
That  is,  assume 

P = (JTJ)"1 

(2.4-12) 


T-' 

E = J J 


Afl>A 

F - J J 


Then , 


T 

J1  J 

Arp 

J1  J 


rn  a 

J J 

Aip  A 

J1  J 


-1 


,~1 


E 

F 


-1 


P ( I+ESE  P ) -PES 
-SETP  S 


(2.4-13) 


where , 

T -1 

S * (F-E  PE) 

A fp  A A W rp  A v 1 

* (J  J - J JPJ  J)  (2.4-14) 

/\ 

The  new  solution  for  B is  found  as  follows  using  eq.  (2.4-11) 
and  eq.  (2.4-13) 

A rp  rp  Aip 

B - -SE  PJ  Z + SJ  Z 

Am  m 

(2.4-15) 


a m t 

SJ  (Z  - JPJ  Z) 


SJ1 (Z  - JA) 


A 

The  solution  A is  found  in  a similar  way,  i.e. 


-r.-- irrrrr.  ■mzi:,  •.  j - 
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Case  3: 

Suppose  the  order  of  the  polynomial  h is  chosen  to 
equal  g,  and  whenever  an  additional  measurement  is  taken  at 
next  maintenance  interval,  the  order  of  the  polynomial  will 
be  increased  by  1.  Then,  eq.  (2.3-4)  becomes: 


TV 

j=o 

kT-1  . 

Z j1. 
j~o 

^.g! 

. . Z j j 

j*0  i 

TV1 

j=o 

a0  ' 

"1_C1  ' 

"TV 

j=0 

2kT-l  . 

Z j1 
j=o  * 

1 

2kT-l  ! 

I j9! 
j-0  | 

i 

i 

2KM.  . 
J 

j=o 

al 

• 

1“C2 

• 

gkT-i  0 

Z 

j*o 

E 3 . 
j=0 

i 

gkT-l  ! 

s j9! 
j=o  | 

i 

• • 

gkT-l  . 

Z j9*1 

j=o 

• 

• 

ag 

• 

l‘c9 

(g+l)kT-l  0 

E j 
j=° 

(g+1) kT-1  . 

Z j1 

j=o  * 

(g+1)  kT-1  j 

1 3 ! 

• • , _ | 
3=0  ! 

1 

(g+1) kT-1  , 

Z j9*1 
j=0 

ag+l 

or  (2.4-17) 


J 

! H ■ 
» 

_ A , 

A 

Z 

. J 

1 

! G J 

A 

B J 

A 

. z . 

where  J,  Z were  defined  in  eq.  (2.3-3),  J,  A,  Z were  defined 
in  eq.  (2.4-2)  with  h = g,  and, 


!? 


kT-1 

_ .< 

S 3 


H = 


2kT-l 
E j' 
j=0 


^ 3 


G = (^+1)f-1jg+l 
j =0 


and  B = a. 


That  is  the  matrices  A and  B can  be  expressed  as 


J j H' 


(2.4-19) 


A similar  algorithm  as  used  in  case  2 can  be  employed  (24, 


25) . That  is, 


J Hi 


J GJ 


where  S * (G  - JJ^H)-1 


j“1(I+HSJj“1)  -J^HSl 


-SJJ 


(2.4-20) 


i. 


W 


CHAPTER  3 


I 

• FAULT  PREDICTION 

For  the  fault  prediction  algorithm  it  is  assumed  that 
the  value  of  a component,  drifting  due  to  shock  damages,  is 
known  at  a fixed  set  of  points  in  the  "time  domain".  Using 
this  data,  one  estimated  the  unknown  failure  dynamics  for  the 
^ parameter.  This  is  then  used  to  compute  the  component's  life- 

time in  the  "shock  domain";  that  is,  the  number  of  shocks 
required  to  cause  a failure. 

^ In  the  following  paragraphs  some  formulas  and  proper- 

ties of  component's  statistics  will  be  developed.  Using  esti- 
mated lifetimes,  some  decision  rules  will  be  introduced  to 
^ compute  the  optimal  time  at  which  to  replace  the  component  so 

as  to  minimize  the  cost  function.  Several  terms  are  defined 
as  follows: 


» 


» 


» 


» 


(1)  L 

(2)  Tr 

(3)  Pf 

(4)  Pr 


The  component's  estimated  lifetime  in  the  shock 
domain;  i.e.  the  number  of  shocks  required  to 
cause  a failure. 

The  optimal  time  at  which  to  replace  the  component 
so  as  to  minimize  the  cost  function. 

The  probability  that  the  failure  occurs  before  the 
replacement  time  Tr« 

The  probability  that  the  component  is  replaced 
by  a new  component  at  replacement  time  Tf;  i.e. 
the  probability  that  the  component's  lifetime  is 
longer  than  Tr« 
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(5)  fL(t): 


(6)  Tf  : 

(7)  T : 

(8)  T*  : 


The  probability  density  function  that  the 
component  receives  the  Lth  shock  at  time  t, 
where  0 < t <.  T^. 

The  expected  lifetime  of  components  which 
failed  before  replacement. 

The  expected  lifetime  of  components  which 
either  failed  before  Tf  or  were  replaced  at  Tr 
The  expected  lifetime  of  components?  i.  e. 


T*  a Tl  T -►  oo 

J r 

Let  p^(t)  represent  the  density  function  of  the  Poisson 
distribution  with  parameter  (kt)  and  E_ (t)  represent  the 

ii 

corresponding  distribution  function;  i.e. 


pi(t)  = ■ e'kt  i-0,  1,  2,  (3-1) 

L-l 

E_  (t)  * Z p. (t)  (3-2) 

L i-0  x 


Several  properties  can  be  stated  as  follows: 


(prop  1)  P = E_  (T  ) 
r Li  r 

(3-3) 

L-i  L 

P_  - E (kt)1  _ -kT_ 

r i-0  "I~l  e r 

L-l 

- 2 Pi(T_) 

i-0  1 r 

-EL(Tr) 

(prop  2) 


1 - EL(Tr) 


(3-4) 
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» S P±(T_) 
i-L  1 r 
L-l 

- 1 - E P,  (T  ) 
i-0  1 r 

» 1 - E_  (T  ) 
ii  r 


(Prop  3)  P + Pf  - 1 


(3  -5) 


(Prop  4)  ^.Tr 
0 


/„r  Pi(t)  at  - i u - Ei+1(Tr>> 

fT]C  Pi<t)  dt  - fr 
J o 1 J n 


•T“  (kt)1 


0 

k1 

IT 


e_kt  dt 


fT]C 

J n 


e"kt  dt 


(3  -6) 


(3  -7) 


Using  the  identity  eq.  (3-8) 


/*” 


ax  , 
e dx 


eax  E 
r-0 


Eq.  (3-7)  becomes 


(-1) 


ml  x 


m-r 


(m-r)  1 a 


r+T 


(3  -8) 


fr 

0 


pi(t)  dt 
i 


-kt 


kJ 

IT 


1 

"k 


t (-l)r  i 1 t 


i-r 


r-0 


, -k.  0 i 1 

{ * 7T+T 

k 


{ 1 - e”kTr  E 


(i-r)  1 (-k) 

-kT  1 11  Tr1_r 

» Kir  E r 


r+1 


(kTr) 


r-0  (i-r) ! k 
i-r 


r+T 


r-0  (i-r) 1 


} 
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i 

L 

j=>0 


(kT  )j 
r 

T ~ 


Z Pi <T_) } 
j-0  3 r 


- -T  {l  - WV1 

(Prop  5)  f.(t)  = (3-9) 

\ (1-el(V» 

where  0 < t < T 
— r 

If  the  interval  (0,  Tr]  is  divided  into  N subinter- 
vals, it  is  seen  that 

Prob  [ ( i— 1) A < t < i A]  * (3-10) 

N 

T 

where  A - 

N 

Assume  A is  small  enough  such  that  the  probability  of  more 
than  one  shock  in  the  subinterval  ((i-l)A,  iA]  is  equal  to 
0.  The  probability  of  one  shock  in  the  subinterval  ((i-l)A, 
iA]  is  equal  to  kA,  where  k is  the  constant  rate  of  Poisson 
shocks. 

For  a given  t in  the  subinterval  ((i-l)A,  iA],  the 
probability  of  the  shock  in  this  subinterval  is: 


» 


Prob  {I,  shock  occurs  in 


((i-l)A,  iA  / ( (i-1) A < t < iA)> 


Prob{(L-l)  shocks  occur  in  (0,  (i-1)  A]} 


* Prob  {one  shock  occurs  in  ((i-1) A,  iA]} 


[k(i-l) 


-k(i-l)A  * 


(L-l)  ! 


pL-l ( (i-1) A)  * kA 


(3-11) 


From  Baylar  theory. 


Prob  {have  Ltn  in  ((i-1) A,  iA]} 


Prob{ (i-1) A<tiiA}* 


I Prob{  ( j-1)  A<ti.jA}* 

j-1 

*Prob{Lt^1  shock  occurs  in  ( (i-1)  A<tjiA/  (i-1)  A<t<iA} 


*Prob{Lt  shock  occurs  in  ( ( j-1) A<t< jA/ ( j-1) A< t< jA} 


I PL-1(I1~1U)  ** 
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PL-l((i_1)A) 


N 

1 PL-l((j“1)A> 

j-1  L1 


(3-12) 


Since  A is  very  small,  two  approximations  can  be  made  as 
follows: 


(i)  (i-l)A-t 

(ii)  Left  hand  side  of  eq.  (3-12)  is  approached  to 


That  is,  eq. 


fL(t) 

(3-12) 


fL(t) 


* A 


can  be  rewritten  as 


PL-i (fc) 


N 

Z PT  t ( (j-1) A) A 

j-i  L 1 


If  A -*■  o,  by  definition  of  the  integral  the  summation 
becomes  an  integral,  such  that 


PL-l(tl 

rp  a /s 

' r pL-1(t)  dt 
0 


(3-14) 
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From  (Prop  4)  , eq.  (3-14)  becomes 


f L (fc)  * 


Pt.- 


L-l(t) 


K {1-el<V} 


(Prop  6) 


T = L * 1 ~ EL-H(Tr} 
if"  K 


1 - E.  (T  ) 
L r 


Since  ^ is  the  expected  lifetime  of  the  components  which 
failed  before  replacement. 


Jf  "(  r t fL(t)  dt 

/n 


(r  4^1.  dt 


t * (kt)L_1  _-kt  _ 

(L-l)  ; G dt 

f (i  - w 


L r r (kt)1'  _ -kt 

Ll e dt 


£ {1-EL(Tr)} 


L * 


f 


PL(t)  dt 


{1  - EL(Tr) } 


(3-16) 
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From  (Prop  4),  eq.  (3-16)  becomes 
1 _ L . C1  " ELtl(Tr’> 

^ £ 1-  ll  m /mV! 


f k 


- E (T) 
L r' 


(Prop  7)  T = — { 1 - E . (T  ) } + T ET  (T  ) 

^ L+l  r 1 r L r 


T * P-  T,  + P T 
f f r r 


(3-17) 


L 1 - E . (T  ) 
{i  - EL(T  )>*  ^ ) 

L r 


* — (1  * EL+l(Tr)}  + Tr  W 


+ T Et  (T  ) 
r l r 


* L 
(-Prop  8)  T = — 
k 


(3-18) 


Ei(Tr)] 


i-1  (kt)j 


= I 


.-kt] 


Tr=  00  j=0  j! 


i.  e.  The  probability  of  having  a finite  number  of 
shocks  in  an  infinite  time  period  is  zero. 


* A 
T = T] 


Tr=°°  k 


— {1  - el+1  <V  )+Tr  VTr,] 


(Prop  9) 


d(P.) 

(1)  — — -PL-l(Tr) 

d(kTr) 


(3-19) 


kTf  - L 

d (kTf ) 
d(kTr) 


'L+l'  r' 


1 - EL(Tr) 


l,£W± 


pLITr>  - {^WV 

{1  - E. (T  ) }2 

Li  2. 


PL-1 


(T 


d (kT) 

(Prop  11)  - E_  (T  ) 

d(kT  ) L r 


T 


L 

— (1 
k 


- E 


L+l 


kT= 


d (kT) 
d(kTr) 


L {1  - EL+1 
L <PL(Tr)} 


(3-22) 


(T  ; } + T E_  (T  ) 

r r L r 

(Tr)}  + kTrEL(Tr) 

+ kTr(-pL_1(Tr) ) + EL(Tr) 


Since, 


- t PL(Tr)  - XVWfV  + EL(Tr) 
<kT  >L  -kT 

L p.  (T  ) = L * £ — e r 

L r LI 


d (kT) 
• # “ 


(kT  ) 


(kT  ) 


L-l 


(L-l)  i 


-kT 
e r 


kTr  PL-1(V 


W 


a(kT. ) 


Once  the  estimated  failure  time  and  statistics  of  the 
component  under  study  have  been  computed  it  remains  to  make 
a replacement  decision.  This  should  be  based  on  the  proba- 
bility of  the  component's  failure  "P^",  the  cost  of  replac- 
ing the  component  "C  ",  and  the  cost  of  an  on-line  failure 
"Cf".  That  is,  the  cost  function  to  be  optimized  can  be 
expressed  as  either 

1 

Cost  - — (CR  + CfPf)  (3-23) 

or 

Cost  = CfPf+  cRpr  (3-24) 

One  possible  replacement  criterion  is  based  on  the 
cost  of  an  on-line  failure  and  average  lifetime  wastage  due 
to  replacing  the  component  before  its  actual  failure. 


Cost  * CfPf  + Cw(kT  -kT) 


(3-25) 


where  Cf  and  Cw  are  two  weight  constants.  The  first  term 
of  eq.  (3-25)  represents  the  cost  of  an  on-line  failure  and 
the  second  term  represents  the  lifetime  wastage.  The  re- 
placement time  Tr  which  simultaneously  minimizes  the  cost 
of  an  on-line  failure  and  lifetime  wastage  should  satisfy 
eq. (3-26)  or  eq. (3-27) . 


“-Cf  Pl-i'V  - cw  W 


(3-26) 


■ n pvp 
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therefore. 


C,  (kT  ) 

r _ i \ * r 


L-l 


(—  - 1)  * 

(L-l) ! 


L-2  (kT  ) 
£ — -£- 

i=0  i! 


(3-27) 


Table  (3-1)  gives  the  replacement  time  kTr  with 


different  weight  constants  (Cf/Cw)  and  different  estimated 


lifetimes.  Fig.  (3-1)  shows  Tr  curves  being  a function  of 


L with  different  (C^/Cw) . 


1 


replacement  time  (kT  ) for  different  (C-/CL.)  and  estimated  lifetimes  (L) 
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CHAPTER  4 
SIMULATIONS 


4.1  Poisson  Shocks  Generator 

To  generate  a series  of  Poisson  Shocks  for  the 
simulation  of  an  on-line  maintenance  process  by  computer, 
uniformly  distributed  random  numbers  are  needed.  A linear 
congruential  method  (26,  27)  can  be  used  to  generate  num- 
bers uniformly  distributed  on  the  interval  (0,  1)  by 
computer.  That  is,  if  the  following  relation  exists  bet- 
ween a set  of  numbers  {x^} 


xn+l  » <*n  *<2  +1>  +b)modm 


(4.1-1) 


where 


and 


s * Q/2,  m = 2U 

Q = the  number  of  shift  register  bits 
in  the  computer 
b = any  odd  number 


then,  the  sequence  (x^  is  of  period  m and  0 <.  x^  < m-1. 

If  each  x^  is  adjusted  to  the  range  (0,  1)  by  dividing 
by  m;  i.e.,  u^  * x^/m,  then  u^  is  a sequence  of  uniformly 
distributed  random  numbers  between  0 and  1.  That  is, 

p0(u)  = 1 (4.1-2) 

where  0 i u « 1 

In  order  to  convert  the  uniform  distribution  U to  a Poisson 
distribution  N,  assume: 


40 

y = - | In  (1-u)  (4.1-3) 

u=l-e-k*  (4*1-4) 

where  y > 0 

and  k is  a positive  constant. 


By  the  change  of  variable  method  (28,  29) , 
density  function,  py(y),  becomes 

PY(y)  - P„(u)  *|  — | 

Y u dy 


- 1 


* ke 


-ky 


the  probability 


(4.1-5) 


Let  y represent  the  time  interval  between  shocks,  i.e., 
the  times  of  successive  shocks  are  y^  , y2»  ....  etc.  Then, 
the  distribution  of  the  number  of  shocks,  n,  in  a period 
T can  be  solved  as  follows: 

Prob(N=n)  = Probfy^yj  +....+yn+1  ^ T) 

(4.1-6) 

-Prob(y1+y2+ +Yn  ^ T) 


From  eq.  (4.1-5)  it  can  be  seen  than  Y is  distributed  as 

12  m 

(X  with  2 degrees  of  freedom)  and  Z y,  is  distributed 
* j»l  j 

12 

as  ^^(x  with  2m  degrees  of  freeedom)  . Therefore,  (30) 
eq.  (4.1-6)  becomes 


Prob (N*n) 


Prob(x22  (n+i)  2kT)  ” Prob(x22n  2 2kT) 
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4.2  Simulation  of  a Perfect  System 

The  computer  simulation  of  a perfect  (no  measurement 
errors,  and  identical  drifting  parameters)  on-line  mainte- 
nance system  can  be  performed  using  the  following  steps: 

t 

Step  1 Define  Constants 

(1)  T:  maintenance  interval 

(2)  k:  time  constant  for  Poisson  Shock  i.e 

average  number  of  shocks  per  unit 
time. 

(3)  xQ:  seed  for  uniform  distributed  random 

number  generator 

(4)  Cf/C^:  the  ratio  of  two  weight  factors 
Step  2 Generate  Poisson  Shock 

(1)  Generate  Poisson  Shock  y^,  y2 , y3.... 

(2)  Calculate  number  of  shocks  received  in 
each  maintenance  interval  n^,  n2 , n3*. 

Step  3 Define  Values  of  Component  Parameter 

(1)  C(0)  = 1 

(2)  C(i)  = a.^  where  i is  an  integer  and 

0 < i < L and  0 < a^  <,  1 

(3)  C(i)  ■ 0 where  i is  an  integer  and 

1 > L — L is  the  lifetime  of  the 


component  in  the  "Shock  Domain". 
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Step  4 Observe  Component  Performance  at  Next 
Maintenance  Time  mT 

(1)  At  maintenance  time  = mT,  the  total 
number  of  shocks  received  are  N 
therefore  the  performance  data 


m 

l n., 
i*l  1 


m 


C (N) 


0,  then  an  on-line 


(2)  Case  Is  If  C 

m 

failure  occurs. 

GO  TO  STEP  9 

Case  2:  If  Cm  > 0,  then  the  component 

performs  properly. 

GO  TO  STEP  5 

Step  5 Calculate  Failure  Dynamics 

The  new  set  of  constants  {a^}  of  the  failure 
dynamics  can  be  calculated  from  eq.  (2.4-23) 
and  eq.  (2.4-24) 

Step  6 Solve  Component ' s Lifetime 

The  lifetime  of  the  component  is  then  solved 
from  eq.  (2.1-5) . 

Step  7 Estimate  Component 1 s Replacement  Time 

The  replacement  time,  Tr,  can  be  predicted 
from  eq.  (3-27) . 

Step  8 Make  a Decision 

Case  1 : If  Tf  < (m+l)T,  then  the  component 

will  not  be  replaced. 
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GO  TO  STEP  4 

Case  2:  If  Tf  > (m+l)T,  then  the  component  will 
be  replaced. 

GO  TO  STEP  9 

Step  9 Replace  component 

Replace  the  component  and  reset  the  count 
m to  zero 

GO  TO  STEP  4 


. j 


t 
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Example : 

The  computer  simulation  of  a on-line  periodic  main- 
tenance system  is  performed  for  600  maintenance  intervals. 
Assume  the  system  is  subjected  to  a Poisson  Shock  with  con- 
stant k»0.1  shocks/hour  and  each  interval  T is  20  hours. 

The  normalized  values  of  a drifting  parameter  due  to  shocks 
deunage  are  as  followings: 


C(0)  - 

1.0 

C (14) 

a 

0.9 

C(l)  - 

0.998 

C (15) 

a 

0.89 

C (2)  - 

0.995 

C (16) 

a 

0.88 

C (3)  - 

0.989 

C (17) 

a 

0.87 

C (4)  - 

0.982 

C (18) 

a 

0.85 

C (5)  - 

0.975 

C ( 19 ) 

a 

0.83 

C (6)  - 

0.968 

C (20) 

a 

0.8 

C (7)  » 

0.96 

C (21) 

a 

0.75 

C (8)  * 

0.952 

C (22) 

a 

0.68 

C (9)  - 

0.944 

C (23) 

a 

0.6 

C(10)« 

0.936 

C (24) 

a 

0.5 

C(ll)« 

0.927 

C (25) 

a 

0.35 

C (12) ■ 

0.918 

C (26) 

a 

0.2 

C (13) ■ 

0.909 

C (27) 

a 

0.05 

and,  C(N)  - 0 if  N >28 


, 


Table  (4.2-1)  gives  total  number  of  replacements  and  failures 
within  600  maintenance  intervals  with  different  values  of 
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C*/CT 


'f  W.  Table  (4.2-2)  shows  the  total  number  of  replacements 
and  failures  with  different  replacement  times. 


Since  the  cost  function  is 


Cost 


Cf*  Pf  + Cw*  (kT  - kT) 


the  overall  cost  can  be  expressed  as 


Cost 


= if  * 


(No.  of  failures) 


'W 


+0.1  * (280* (No.  of  components  used)  - 12000) 


The  overall  cost  with  different  methods  and  different  values 
of  C^/Cw  are  given  in  Table  (4.2-3). 
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Table  (4.2-1) 

Total  replacements  and  failures  within  600 
maintenance  intervals  for  different  Cf/Cw 


! 


Table  (4.2-2) 

Total  replacements  and  failures  within  600  maintenance 
intervals  for  methods  of  constant  time  replacement 


Constant 
replacement  time 

No.  of  replacement 

No.  of  failure 

every  6 

intervals 

100 

0 

every  7 

intervals 

85 

0 

every  8 

intervals 

75 

0 

every  9 

intervals 

65 

1 

every  10 

intervals 

59 

1 

every  11 

intervals 

48 

6 

every  12 

intervals 

39 

11 

H 


Table  (4.2-3) 

Overall  cost  with  different  methods  and  different  ®f/Cw 


v <vtw 

Costas. 

50 

75 

100 

150 

200 

Methods 

every  6 
intervals 

1600 

1600 

1600 

1600 

1600 

every  7 
intervals 

1096 

1096 

1096 

1096 

1096 

every  8 
intervals 

900 

900 

900 

900 

900 

every  9 
intervals 

698 

723 

748 

798 

848 

every  10 
intervals 

530 

555 

580 

630 

680 

every  11 
intervals 

612 

762 

912 

1212 

1512 

every  12 
intervals 

750 

1025 

1300 

1850 

2400 

the  algorithm 

690 

471 

512 

668 

763 
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4.3  Simulation  of  an  Imperfect  System 

The  computer  simulation  of  an  imperfect  on-line  main- 
tenance system  can  be  performed  by  modifying  the  previous 
algorithm.  Instead  of  using  the  seine  set  of  C(i)  for  every 
component,  a noise  is  added  into  C(i)  such  that  the  lifetime 
of  each  component  may  not  be  the  same  but  the  average  life- 
time is  still  unchanged. 

The  total  number  of  replacements  and  failures  with 
= 100  in  600  maintenance  intervals  are  solved  by  using 
different  methods.  The  results  with  different  noise  levels 
are  showed  in  Table  (4.3-1). 

Table  (4.3-2)  gives  the  overall  cost  with  different 
methods  and  different  noise  levels. 
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Table  (4.3-2) 

Overall  cost  for  different  methods  at  different 
noise  levels 


«s^\|23 

0 « 

20  « 

30  % 

40  % 

60  % 

method 

every  6 
intervals 

1600 

1600 

1600 

1600 

2200 

every  7 
intervals 

1096 

1096 

1096 

1280 

2008 

every  8 
intervals 

900 

900 

1200 

1300 

2128 

every  9 
intervals 

748 

848 

948 

1376 

2432 

every  10 
intervals 

580 

880 

1380 

1980 

2364 

every  11 
intervals 

912 

1340 

1340 

1340 

2452 

every  12 
intervals 

1300 

1728 

1828 

1984 

2612 

the  algorithm 

512 

752 

980 

752 

1608 
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CHAPTER  5 
CONCLUSION 

In  the  proceeding  chapters  we  have  described  a fault 
prediction  algorithm  for  on-line  maintenance  systems.  This 
algorithm  essentially  can  be  separated  into  three  major  steps: 
1)  applying  a curve  fitting  technique  to  solve  the  failure 
dynamics  by  using  past  performance  data,  2)  predicting  the 
lifetime  of  the  system  in  the  shock  domain  from  the  failure 
dynamics,  3)  estimating  the  replacement  time  which  simultan- 
eously minimizes  the  cost  of  an  on-line  failure  and  lifetime 
wastage.  In  addition,  sequential  refinement  schemes  have 
been  developed  to  solve  the  problem  of  inverting  a potentially 
high  dimensional  matrix.  Thus  without  having  to  repeat  the 
entire  calculation,  the  new  sets  of  failure  dynamics  can  be 
obtained  recursively  based  on  the  old  estimates  and  new  data. 

The  algorithm  has  been  tested  in  a variety  of  situa- 
tions such  as:  perfect  and  imperfect  system,  different  levels 
of  noise,  different  sets  of  Poisson  shocks.  The  results  have 
been  found  to  be  surprisingly  effective  in  predicting  failures 
with  relatively  little  wastage  of  lifetime  and  on-line  fail- 
ure cost. 

Finally,  similar  algorithms  for  the  replacement  cri- 
terion based  on  the  cost  functions  of  eq.  (4.1-23)  and  eq. 
(4.1-24)  have  been  studied  using  the  properties  introduced 
in  chapter  3.  These  algorithms  yielded  very  good  results. 
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//  FOR  MAX  01  JUN  77  23.409  HRS 
♦LIST  ALL 

♦IOCS  (CARD, DISK,  1443  PRINTER) 

♦CNEWCRD  INTEGERS 
♦NCNPROCESS  PROGRAM 

SUBROUTINE  MAX(IOODE,AM,BM,C34,M,N,K) 

C. 

C.  THIS  IS  THE  SUBROUTINE  FOR  MATRIX  OPERATIC!! 
C.  AM(M,N)  EM(N,K)  CM(M,K) 

C.  IOOCE-1  CM=AMffiM 

C.  IOCCE-2  QWM-BH 

C.  IOODE-3  OKAM*EM 

C. 

DMNSICN  AM (20, 20)  ,BM(20,20)  ,CM(20,20) 
IF (IOOOE-2  ) 100,200,300 
iuO  CONTINUE 

DO  10  I»1,N 
DO  10  J»1,N 

(34(1, J)  =*  AM(I,J)  + BM(I, J) 

10  CONTINUE 

GO  TO  400 
200  CCNTINUE 

DO  20  1*1, N 
DO  20  J-1,N 

CM(I,J)  =*  AM(I,J)  - HM(I, J) 

20  CCNTINUE 

GO  TO  400 
300  CONTINUE 

DO  30  I»1,M 
DO  30  J*1,K 
CM(I,J)  * 0.0 
DO  5 Ir=l,N 

04(1, J)  - (34(1, J)  + AM(I,L)*BM(L,J) 

5 CCNTINUE 

30  CONTINUE 

400  CCNTINUE 

RETURN 
END 

//  FOR  FNOT  01  JUN  77  23.412  HRS 
♦IOCS  (CARD, DISK,  1443  PRINTER) 

♦LIST  AIL 

♦NCNPROCESS  PROGRAM 
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— 


; 


I 


DIMENSION  NSK(600) , IDESN(600) , QB(20) 

EQUIVALENCE  (FIN  (1,1)  ,TSHUK(1)  ) 

E)QUIVAIENCE(FINEW(1,1)  ,TSHUK(401)  ,DUM2(1,1) ) 

EQUIVALENCE  ( DUM1(1,1)  ,TSHUK(801)  ) 

EQUIVALENCE  (DUM(1,1)  ,TSHUK(1201)  ) 

EQUIVALENCE  (A (1,1)  ,TSHUK(1601) ) 

EQUIVALENCE  (H(  1,1)  ,TSHUK(1621) ) 

EQUIVALENCE  ( E (1,1) ,TSHUK (1641) ) 

EQUIVALENCE  (C  ( 1, 1)  ,TSHUK(  1661)  ) 

EQUIVALENCE (OB ( 1) ,TSHUK(1681) ) 

EQUIVALENCE (CM(1,1)  ,TSHUK(1701) ) 

EQUIVALENCE (AM(1,1)  ,TSHUK(2101) ) 

EQUIVALENCE (EM(1,1)  ,TSHUK(2501) ) 

EQUIVALENCE  (CH(1  ) ,TSHUK(2901) ) 

CR=5 

PR=6 

C. 

C.  THIS  SECTION  IS  TO  (3NERA1ED  POISSON  SHOCK 
C. 

C.  SEED  THE  SEED  FOR  THE  R.  N.  GENERATOR 
READ (CR, 9101)  SEED 
9101  FORMAT (F  5.0) 

C.  NFND  NLM  OF  RANDOM  NUMBER  TO  BE  GENERATED 
READ (CR, 9101)  NFND 
9101  PQRM\T(14) 

C. 

C.  INPUT  THE  OBSERVATION  PERIOD  TPD 
C. 

READ  (CR,  901)  TPD 
901  FORMAT  (F15. 5) 

C.  TCNST  TIME  CONSTANT 

READ  (CR,9101)  TCNST 
CNST=TPD/TCNST 
WRITE  (PR,  9106)  SEED 

9106  FORMAT (1H  , 'THE  SPm  FOR  THE  RANDOM  NO.  IS  ' ,F5.0  ) 

WRITE  (PR,  9107)  NFND 

9107  FORMAT (1H  , 'THE  NO.  OF  POISSON  SHOCK  TO  BE  GENERATED  IS 
WRITE  (PR,  9108)  TCNST 

9108  FORMAT (1H  ,'AVE.  NO.  OF  SHOCK  PER  OBSERVATION  PERIOD  IS 
C. 

C.  X (N+l)  = (X  (N)  *AN  +EN)  MOD  CN 
C.  X(N)  IS  R.N.  N(0, 16383) 

C. 

DO  9115  INRNEKL,NFND 

AN-129.0 

BN*  111.0 

CN-16384.0 

SI-SEED* AN  +BN 

IS-S1/CN 

S2»S1-IS*CN 


t 
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SEEIX>2 

! § c. 

C.  SEER  NORMALIZED  R.  N.  N(0,1) 
C.  TSHUK  IS  POISSCN  SHUCK 
C. 


SEEte»SEED/CN 

TSHUK (INFND)  » (-  CNST)  *ALOG(  1-SEER) 


946  PCW«r(lH  ,10  (15, 2X)) 

950  CONTINUE 
ICNP-5 

DO  100  101,12 
ICNP-ICNP  +1 
WRITE  (PR,  101)  ICNP 
101  FORMAT  (/////, ' PERIOD  ' ,14) 

ICNPT-0 

ISUM-0 

DO  105  I-1,MAXPT 
ICNPT-ICNPT  +1 
ISl*MSUMfNSK(I) 

IF (ISOM-28)  110,120,120 

no  if  (icnpt-icnp)  105,  m.  Ill 

in  WRITE  (PR,  115)  I 

115  FORMAT ('  AT  TEST  POINT  ',  I 5,  'REPLACE') 

ICNPT-0 

ISURH) 

00  TO  105 

120  WRITE  (PR,  125)  I 

125  FORMAT ('  AT  TEST  POINT',  15,'  FAIL') 

ISOM-O 

ICNPT-0 

105  CONTINUE 
100  CONTINUE 
C. 

C.  INPUT  THE  CHARACTERISTIC  OF  THE  DEVICE 
C. 

C,  C(Q)  » 1 

C.  C(N+1)  ALWAYS  LESS  THAN  C(N) 

C.  4 DATA  CARDS  WITH  10(F6.5  IX) 

WRITE  (PR,  9315) 

9315  FORMAT  (1H  , "THE  CHARACTERISTIC  OF  THE  DEVICE  ARE  ') 
READ (CR, 9310)  I 
9310  FORMAT  (78X,A1) 

DO  9300  1-1,4 
IL>-(I-1)*10  +1 

iH-mo 

READ  (CR,  9350)  (CH(J)  ,J-LL,LH) 

9350  FORMAT  ( 10(F6.5,1X)) 

WRITE  (PR,  9360)  (CH(J)  ,J-LL,LH) 

9360  FORMAT (1H  ,10(F7.5,1X) ) 

9300  CONTINUE 
C. 

C.  INPUT  THE  RATIO  OF  CF/CW 
C.  CODE  FOR  IDESN  ARE 
C.  IDESN (IDE)  - 1 FAIL 

C.  IDESN  (IDE)  - 2 CONTINUE 

C.  IDESN (IDE)  * 3 REPLACE 


c. 

MT-TCNST 

READ  (CR,  3002)  CFW 

3002  FORMAT  (F15. 5) 

WRITE  (PR,  3003)  CFW 

3003  F0RMAT(1H  , 'THE  RATIO  CF/CW  IS  ' ,F15.5) 

MT-TCNST 

IEE=0 

3005  ICT»0 

3006  IDB=IDE+1 

IF  (IEE-MAXPT)  3010,3010,4000 
3010  IC3MCIM-1 

IF(ICT-l)  3015,3015,3115 
3015  ITSK”  NSK(IDE) 

IF(CH(rrSK) ) 3020,3020,3025 

C.  

C.  FIRST  TIME  TEST  AND  FAIL 
C. 

3020  WRITE  (PR,  3022)  IEE.ITSK 
3022  FORMAT (1H  ,"AT  TEST  POINT'  ,I5,3X, 
l'THE  DEVICE  TOTAL  RECEIVED'  ,15, 2X, 

2 'SHOCKS  AND  FAIL  TO  OPERATE') 

H3ESN(IEE)*1 
GO  TO  3005 
C. 

C.  FIRST  TIME  BUT  NOT  FAIL 
C. 

3025  F2N(1,1)=1.0/MT 

CB(ICT)  = CHdTSK) 

C(l,l)  » 1.0  -OB (Id) 

NORED=l 

A(1,1)=FIN(1,1)*C(1,1) 

GO  TO  502 
C. 

C. 

C.  THIS  IS  THE  ROUTINE  TO  SOLVE  THE  PARAMETERS 

C NORED  THE  ORDER  OF  THE  COEFFICIENTS 

C.  FI  THE  INVERSE  CF  THE  MATRIX  F 

C.  F X A » C 

C.  EH  * A « C 

C.  EG  * B « D 

C.  H (NORED  , 1) 

C.  E (1  , NORED) 

C.  G , D IS  CONSTANT 
3115  CONTINUE 

ITSK=nSK  + NSK(IEE) 

IF  (CH(ITSK) ) 3020,3020,1000 
1000  CONTINUE 
C. 

C.  CALCULATE  H C 


■■Hp 


NMEl-NOBEDH 

SlM>0.0 

00  1010  I»1,NQRED 
Ufl^(I-l)  * MT 
IF(IML)  1002,1002,1003 

1002  IML-1 

1003  U«-  I*MT  - 1 
DO  1005  J-  LML  , 

SIM  « SUM  + J **  NORE1 
1005  CONTINUE 

H(I,l)*StM 
C(I,1)-1.0-  OB (I) 

1010  CONTINUE 
C. 

C.  CALCULATE  G 
C. 

IMD*  I^H+1 
LMR=  IML*  MT 
DO  1110  J ■ IML  ,IMH 
SUM*  SUM  + J **  N0RE1 
1110  CONTINUE 
G-StM 
C. 

C.  CALCULATE  E 
C. 

DO  1111  1-2, NOSED 
lin  E(1,I)  » 0.0 

E(l,l)  * (NQRE1*MT) 

U«=  (NOEE1  ) * MT  - 1 
DO  1120  J=  2,NOHED 
DO  U15  1=  1,  I*H 
E(l,  J)  = E(1,J)  + I **  J 
U15  CONTINUE 
1120  CONTINUE 

CB(NDREl)  » CH(ITSK) 

LHl.O-CB(NOREl) 

CAUL  MAX  (3,E,FIN,DUM,1,NQRED,NQRED) 

CALL  MAX  (3,DUM,H,DIM1,1,NQRED,1) 
S-1.0/(G-D0M1(1,1)) 

FINEW(N0RE1,N0FE1)  =S 
DO  1500  J«1,NQFED 
FINEM(N0RE1,J)  - -S  *D0M(1,J) 

1500  CONTINUE 

DO  1510  J-l,NORED 
DUM(1,J)  - DtM(l,J)  *S 
1510  CONTINUE 

CALL  MAX  ( 3,FIN,H,DUM1,N0RED,NDRHD,1) 

CAUL  MAX(3,DUMl,DUM,DUM2,NC®ED, l*NOFED) 

GALL  MfiX  (1, FIN, DUM2,FINEW,NORED, NOREX), NDRED) 
CALL  MAX  (3,FIN,H,DUM,NQRED,NOKED,l) 


I/.,,., 


DO  1520  I=1,N3REE 
FINEW(I,NQRE1)  « -DUM(I,1)  *S 
1520  CCNTCNUE 

C(NQRE1,1)  *1. 0-QB (NOREl) 

CALL  MAX ( 3 , FINES*/  ,C,A,NOREl,NOREl,l) 

NOHE3>NQFEl 
DO  400  I=1,N0RED 
DO  400  J-1,NQRED 
FIN(I,  J)  » FINEW(I,  J) 

400  CONTINUE 
C. 

C.  FHCM  THE  SET  CP  CCEFFIdENTS  A(I)  SOLVE  THE  LIFE  TIME 
C. 

502  SLMO.O 

WRITE  (PR,  7918)  (A(IT,1)  , HONORED) 

7918  POFMATC  *****  A **'  ,6(F9.4,1X) ) 

J=0 

SU«UMfA(l  ,1) 

IF  (SUM-1.0)  505,510,510 

501  SUM-SLM+A(1,1) 

IF(  NQRED-1)  500,500,513 
513  DO  500  I»2,NORED 
RJ=J 

SDM1=RJ**(I-1) 

SUM«SUMfA  (1,1) * J** (1-1) 

500  OONITNUE 

IF  (S0M-1. 0)  505,510,510 
505  J=d+1 

IF ( J-60 ) 501,501,510 
510  UFE=C+1 

LTFO^UFE-ITSK 
TLIFE*LIFEN/T'CNST 
WRITE  (PR,  7900)  LIFEN 

7900  POMftTC  **  **  **  LIFE=  ’,15,  'SHOCK') 

C. 

C.  USING  NEWTON'S  METHOD  TO  SOLVE  TR 
C.  CFW  CF/CW 
C. 

IF (LIFE  -2)  550,560,565 

550  TR-0.0 

GO  TO  700 

560  Hfr-1.0/(CfW-1.0) 

GO  TO  700 

565  IF  (LIFE-3)  570,570,580 

570  TT^d.O  + SQFT(2.0*CFW-1))/(CEW-1) 

GO  TO  700 

580  IF (LIFE-60) 590 , 567 , 567 
567  TR-9999.0 

GO  TO  700 
590  TR-LIFE*2.0 


— 
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ICQDP=0 
T.TP1  aT.TFR— 1 

UF2=LIFE-2 

FP-1.0 

594  FAT-1.0 

StM=1.0 
TEFto*1.0 
DO  591  1=1,  UF1 
TEPM=TEFM*TRA 
SUM-SUM+TERM 

591  CONTINUE 
F-CFW^IEFM-SUM 

IF (F-0. 0)592,700,593 
593  TR»TR-1.0 

FP-F 
IOODP-1 
GO  TO  594 

592  TR=TR-(F/(FP-F)) 

LIF1=LIFE-1 

UF2-LIFE-2 

610  FAT-1.0 

SUM-1.0 
TEFM-1.0 
DO  600  1=1,  UF2 
TEPM=TEFM*TR/I 
SUM-SUMWERM 
600  CCNTINUE 

IF-CEW*TEPM-SUM 

TER»=TERM*TR/LIF1 

SUM-SUMtTEBM 

F=CFW*TEFM-SLM 

TFN-TR-F/DF 

IF  (ABS(TIW-TR)  -0.001)  700,700,690 

690  TR-TFN 

GO  TO  610 
700  CONTINUE 

TIWR*TPD/TCNST 

'nWR-TTD*ICT 

WRITE  (PR,  710)  HE,  ITSK,OB  (ICT) 

710  FQRAT(£H  ,/,  ’ AT  TEST  POINT'  ,15,  2X, 
l'THE  DEVICE  TOTAL  RECEIVED ' ,I5,2X, 

2' SHOCKS  AND  THE  MEASUREMENT  IS',F15.5) 

WRITE  (PR,  715)  TLIFE,  TR 

715  FORMAT ('  THE  ESTIMATED  LIFE  TIME  IS  ’ ,F9.3,' 
1 ' F7.2) 

IF(TR-TPD)  720,720,730 
C. 

C.  HEPIACE  DEVICE 
C. 

720  NRnE(PR,  722) 


THE  REPLACE  TIME  IS 


FORMAT  (1H+,T7  5,  'IA3ER.  WE  WILL  FEPIACE  IT  NOW') 

U*SN(IEE)-3 

QO  TO  3005 

ITINUE  USING  THE  DEVICE 
WRITE  (PR,  732) 

FORMAT (1H+,T75,  'LA3ER.  WE  WILL  OONTT.  USE  IT') 

H*SN(IDE)»2 

GO  TO  3006 

CONTINUE 


